<!DOCTYPE html>
<html lang="en">
  <head>
    <meta charset="utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1">
    <title>不同分布的二元混合分布 | 阿力木·达依木|Alimu Dayimu</title>
    <link rel="stylesheet" href="https://alim.gitee.io/css/style.css" />
    <link rel="stylesheet" href="https://alim.gitee.io/css/fonts.css" />
    
<link href="//cdn.jsdelivr.net/gh/highlightjs/cdn-release@11.6.0/build/styles/github.min.css" rel="stylesheet">


    


<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/katex.min.css" integrity="sha384-Xi8rHCmBmhbuyyhbI88391ZKP2dmfnOl4rT9ZfRI7mLTdk1wblIUnrIq35nqwEvC" crossorigin="anonymous">
<script defer src="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/katex.min.js" integrity="sha384-X/XCfMm41VSsqRNQgDerQczD69XqmjOOOwYQvr/uuC+j4OPoNhVgjdGFwhvN02Ja" crossorigin="anonymous"></script>
<script defer src="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/contrib/auto-render.min.js" integrity="sha384-+XBljXPPiv+OzfbB3cVmLHf4hdUFHlWNZN5spNQ7rmHTXpd7WvJum6fIACpNNfIR" crossorigin="anonymous"></script>
<script>
    document.addEventListener("DOMContentLoaded", function() {
        renderMathInElement(document.body, {
          
          
          delimiters: [
          {left: "$$", right: "$$", display: true},
          {left: "$", right: "$", display: false},
          {left: "\\(", right: "\\)", display: false},
          {left: "\\[", right: "\\]", display: true},
          {left: "\\begin{equation}", right: "\\end{equation}", display: true},
          {left: "\\begin{align}", right: "\\end{align}", display: true},
          {left: "\\begin{alignat}", right: "\\end{alignat}", display: true},
          {left: "\\begin{gather}", right: "\\end{gather}", display: true},
          {left: "\\begin{CD}", right: "\\end{CD}", display: true},
          {left: "\\[", right: "\\]", display: true}
          ],
          
          throwOnError : false
        });
    });
</script>


  </head>

  <body>
    <nav>
    <ul class="menu">
      
      <li><a href="https://alim.gitee.io/">Home</a></li>
      
      <li><a href="https://alim.gitee.io/post/">Post</a></li>
      
      <li><a href="https://alim.gitee.io/talks/">Talks</a></li>
      
      <li><a href="https://alim.gitee.io/about/">About</a></li>
      
    </ul>
    <hr/>
    </nav>

<div class="article-meta">
<h1 class="title">不同分布的二元混合分布</h1>
<div class="meta">
<span class="author">Alimu Dayimu</span>
<span class="date middot">2022/07/19</span>

<span class="reading-time middot"> 3 min read </span>
<div class="terms">
  
  
  
  
  
</div>
</div>
</div>




<main>



<p>一般情况下，构建模型的时候我们会考虑数据服从单一分布，比如高斯分布、贝塔分布和伽马分布等等。但是有些时候我们会考虑数据可能存在多个分布，比如两个高斯分布，即均值不同，数据出现两个峰值（高斯混合分布），或者是两个不同的分布叠加在一起，比如高斯分布和均匀分布。这种一般可以将不同的分布分别进行拟合，某种程度上也可以认为是一种分类过程，将数据按照分布的不同分成不同的组。在机器学习上比较有名的就是Guassian mixture model，其原理也基本一致。</p>
<p>对于一由有限的<span class="math inline">\(K\)</span>个不同的混合分布组成的分布<span class="math inline">\(f\)</span>表示为
<span class="math display">\[f(x)=\Sigma_k^K\pi_kf_k(x)\]</span>
其混合权重表示为
<span class="math display">\[\Sigma_k \pi_k=1, \pi_k&gt;0\]</span></p>
<p>有了上面的公式以后我们可以通过下面的方法得到一个混合分布的数据：</p>
<ol style="list-style-type: decimal">
<li>以概率<span class="math inline">\(\pi_k\)</span>选取一个分布<span class="math inline">\(f_k\)</span>。</li>
<li>根据上述分布函数生成观察数据。</li>
</ol>
<p>虽然上面是数据生成过程，但是如何确定不同分布的权重<span class="math inline">\(\pi_k\)</span>的同时，估计分布的参数是比较困难的。但是既然我们知道了数据生成机制，那么可以根据这个数据生成机制来构建我们的模型。</p>
<p>接下来部分我们将通过模拟数据并分析来了解这个过程是如何实现的。</p>
<ol style="list-style-type: decimal">
<li>模拟数据，生成由正态分布和均匀分布的混合数据。</li>
<li>用Jags构建分析模型并分析</li>
</ol>
<div id="生成数据" class="section level2">
<h2>生成数据</h2>
<pre class="r"><code>set.seed(2022)
N &lt;- 100
alpha &lt;- 0.3
mu &lt;- 5
sd &lt;- 2
min &lt;- mu
max &lt;- 50
unif_dt &lt;- runif(N, min, max)
norm_dt &lt;- rnorm(N, mu, sd)
latent_class &lt;- rbinom(N, 1, alpha)
Y &lt;- ifelse(latent_class, unif_dt, norm_dt)

hist(Y, breaks=100)</code></pre>
<p><img src="https://alim.gitee.io/post/2022/two-component-mixture-model/index_files/figure-html/unnamed-chunk-1-1.png" width="672" /></p>
<p>上面我们模拟了样本量为100的数据，30%是均匀分布，其余70%是正态分布，其观察到的数据如图所示。</p>
</div>
<div id="模型构建" class="section level2">
<h2>模型构建</h2>
<p>下面是数据分析的Jags模型：</p>
<pre class="r"><code>model &lt;- &quot;
model{
  for(i in 1:N){
    
    # Log density for the normal part:
    ld_comp[i, 1] &lt;- logdensity.norm(Y[i], mu, tau)
    # Log density for the uniform part:
    ld_comp[i, 2] &lt;- logdensity.unif(Y[i], lower, upper)
    
    # The latent class part using dcat:
    component_chosen[i] ~ dcat(probs)
    
    # Select one of these two densities and normalise with a Constant:
    density[i] &lt;- exp(ld_comp[i, component_chosen[i]] - Constant)
    
    # Generate a likelihood for the MCMC sampler:
    Ones[i] ~ dbern(density[i])
    
  }
  
  # Priors:
  probs ~ ddirch(c(1,1))
  lower ~ dnorm(0, 10^-6)
  upper ~ dnorm(0, 10^-6)
  mu ~ dnorm(0, 10^-6)
  tau ~ dgamma(0.01, 0.01)
}&quot;</code></pre>
<p>上面的代码包括两部分，一个正态分布和均匀分布，最后属于哪个分布由伯努利分布来确定。最后虽然我们<code>density</code>那边就结束了，但是为了能够得到正确的分布，我们让<code>density</code>服从伯努利分布，即我们需要最大化后验分布，也就是最大化似然函数。</p>
<pre class="r"><code>dat &lt;- list(N = N,
            Y = Y,
            Ones = rep(1,N),
            Constant = 10)
inits &lt;- list(lower = min(Y)-10,
              upper = max(Y)+10,
              mu = 0,
              tau = 0.01)
jags &lt;- jags.model(textConnection(model), 
                   data = dat,
                   inits = inits,
                   n.chains = 4)
update(jags, n.iter=1000)
samps &lt;- coda.samples(jags, c(&quot;lower&quot;, &quot;upper&quot;, &quot;probs&quot;, &quot;mu&quot;, &quot;tau&quot;), n.iter = 10000 )</code></pre>
<p>上面拟合的模型里面的Constant不能太大，防止density大于1，也不能太小。下面是模型的参数估计，其中<code>tau</code>是精度，即方差的倒数。可以发现下面的参数估计结果跟之前的模拟的结果基本一致。</p>
<pre class="r"><code>summary(samps)</code></pre>
<pre><code>## 
## Iterations = 2001:12000
## Thinning interval = 1 
## Number of chains = 4 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##             Mean      SD  Naive SE Time-series SE
## lower     6.0088 5.59742 0.0279871      0.2865510
## mu        5.5658 0.28707 0.0014354      0.0027789
## probs[1]  0.7271 0.05993 0.0002997      0.0017312
## probs[2]  0.2729 0.05993 0.0002997      0.0017312
## tau       0.2108 0.04327 0.0002164      0.0007612
## upper    44.8976 1.60118 0.0080059      0.0206484
## 
## 2. Quantiles for each variable:
## 
##             2.5%     25%     50%     75%   97.5%
## lower    -2.6416  0.9377  5.3922 11.1039 15.2872
## mu        5.0095  5.3727  5.5620  5.7570  6.1344
## probs[1]  0.5987  0.6886  0.7315  0.7703  0.8311
## probs[2]  0.1689  0.2297  0.2685  0.3114  0.4013
## tau       0.1387  0.1803  0.2067  0.2363  0.3063
## upper    43.4070 43.7943 44.3940 45.4583 49.2235</code></pre>
</div>
<div id="两个混合高斯分布模型" class="section level2">
<h2>两个混合高斯分布模型</h2>
<p>下面的代码是来自<a href="http://doingbayesiandataanalysis.blogspot.com/2012/06/mixture-of-normal-distributions.html">http://doingbayesiandataanalysis.blogspot.com/2012/06/mixture-of-normal-distributions.html</a></p>
<pre class="r"><code>modelstring &lt;- &quot;
model {
    # Likelihood:
    for( i in 1 : N ) {
      y[i] ~ dnorm( mu[i] , tau[i] ) 
      mu[i] &lt;- muOfClust[ clust[i] ]
      tau[i] &lt;- tauOfClust[clust[i]]
      clust[i] ~ dcat( pClust[1:Nclust] )
    }
    # Prior:
    for ( clustIdx in 1: Nclust ) {
    
      muOfClust[clustIdx] ~ dnorm( 0 , 1.0E-10 )
      tauOfClust[clustIdx] ~ dgamma( 0.01 , 0.01 )
      
    }
    pClust[1:Nclust] ~ ddirch( onesRepNclust )
}
&quot;

# Generate random data from known parameter values:
set.seed(47405)
trueM1 = 100
N1 = 200
trueM2 = 145 # 145 for first example below; 130 for second example
N2 = 200
trueSD = 15
effsz = abs( trueM2 - trueM1 ) / trueSD
y1 = rnorm( N1 ) 
y1 = (y1-mean(y1))/sd(y1) * trueSD + trueM1
y2 = rnorm( N2 ) 
y2 = (y2-mean(y2))/sd(y2) * trueSD + trueM2
y = c( y1 , y2 ) 
N = length(y)

# Must have at least one data point with fixed assignment 
# to each cluster, otherwise some clusters will end up empty:
Nclust = 2
clust = rep(NA,N) 
clust[which.min(y)]=1 # smallest value assigned to cluster 1
clust[which.max(y)]=2 # highest value assigned to cluster 2 
dataList = list(
    y = y ,
    N = N ,
    Nclust = Nclust ,
    clust = clust ,
    onesRepNclust = rep(1,Nclust) 
)

model = jags.model(textConnection(modelstring), data=dataList)

update(model,n.iter=1000)

output=coda.samples(model=model,
        variable.names=c(&quot;muOfClust&quot;,&quot;tau&quot;,&quot;clust&quot;,&quot;pClust&quot; ),
        n.iter=1000,thin=1)</code></pre>
</div>
<div id="后记" class="section level2">
<h2>后记</h2>
<ul>
<li>对数据进行描述很重要，学了统计推断以后眼里只有p值，忘记了描述性分析的重要性，分析数据前一定要对数据做描述。</li>
<li>不能简单地直接使用自己掌握的模型，了解数据生成机制很重要，最终的目的是要构建符合数据生成机制的模型。</li>
</ul>
</div>

</main>


  <footer>
  <script src="//cdn.jsdelivr.net/gh/highlightjs/cdn-release@11.6.0/build/highlight.min.js"></script>
<script src="//cdn.jsdelivr.net/gh/highlightjs/cdn-release@11.6.0/build/languages/r.min.js"></script>

<script>
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
</script>

  
  <hr/>
  © <a href="https://alim.gitee.io/">Alimu Dayimu</a> 2021 - 2022 | <a href="https://creativecommons.org/licenses/by-nc-nd/4.0/">BY-NC-ND 4.0</a> | <a href="https://github.com/adayim">GitHub</a> | <a href="https://gitee.com/alim">gitee</a> | <a href="https://alim.gitee.io/index.xml">RSS</a>
  
  
  </footer>
  </body>
</html>
